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Abstract 

Sexually reproducing populations with small number of individuals may go extinct by stochastic 
fluctuations in sex determination, causing all their members to become male or female in a gener- 
ation. In this work we calculate the time to extinction of isolated populations with fixed number 
N of individuals that are updated according to the Moran birth and death process. At each time 
step, one individual is randomly selected and replaced by its offspring resulting from mating with 
another individual of opposite sex; the offspring can be male or female with equal probability. A 
set of N time steps is called a generation, the average time it takes for the entire population to be 
replaced. The number k of females fluctuates in time, similarly to a random walk, and extinction, 
which is the only asymptotic possibility, occurs when k = or k = N. We show that it takes only 
one generation for an arbitrary initial distribution of males and females to approach the binomial 
distribution. This distribution, however, is unstable and the population eventually goes extinct in 
2 N /N generations. We also discuss the robustness of these results against bias in the determina- 
tion of the sex of the offspring, a characteristic promoted by infection by the bacteria Wolbachia 
in some arthropod species or by temperature in reptiles. 



I. INTRODUCTION 



Most species in our planet have small numbers of individuals [lj . Even the human pop- 
ulation, now with more than seven billion people, has gone through periods of very low 
abundances not too long ago [2]. In fact,typical abundance distributions of several groups 
of species has been shown that follow a universal lognormal curve with an excess of rare 
species [U [3] . 

Small communities are prone to extinction for a large number of reasons, such as inability 
to protect themselves or difficulty in finding mates [4|, accumulation of deleterious mutations 
[5] and shear stochastic fluctuations in the environment [&H8] or in the number of males and 
females in the group [9] . If male and female offspring are equally likely to occur, the ratio 
between males and females fluctuates over generations [10] and it might occur that they 
all become male or female, driving the population to extinction. This works against the 
persistence of small communities and it is somewhat puzzling how so many such populations 
do exist. In some species, the females developed ways to control the sex ratio of their progeny, 
producing more females in a male rich environment and vice-verse. This mechanism, termed 
local mate competition JTTJ [12] , may have evolved to avoids extinction by fluctuations in sex 
determination. 

In this work we consider the changes in the sex ratio of a population with fixed number 
of individuals in which a single member is randomly selected to reproduce at each time 
step, being replaced by a male or a female offspring with equal probability. The model 
is a great simplification of the dynamics of real populations, but it captures the random 
character of the process and sheds light on why low abundance species are actually much 
more robust than one could naively expect. In particular, we do not take into account the 
fact that females usually have many offspring and, even if the population size is held fixed, 
there is competition among the offspring and it is the best fit individual who survives. In a 
population consisting mostly of males, it is expected that a female offspring will fare better 
than a male, contributing to balance the sex ratio. We shall not take competition, natural 
selection, spatial structures [13] or aging into account [9], restricting our work to neutral 
evolution, which has been shown to describe the observed universal patterns of abundance 
and diversity, both for hermaphroditic and sexual populations [HHT6]. 
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II. THE DYNAMICAL SYSTEM 



We consider a population with N individuals divided into k females and N — k males. 
We call Pt{k) be the probability of finding k females at time t. The population is updated at 
discrete time steps similarly to the birth and death process proposed by Moran in population 
genetics [17J: at each time step a random individual is selected to reproduce with an available 
member of the opposite sex; after reproduction the selected individual dies and is replaced 
by the offspring, which can be male or female with equal probability The states -P(O) and 
P{N) are absorbing states corresponding to all males and all females respectively and mark 
the extinction of the population, since reproduction becomes impossible when all members 
have the same sex. In a state with k females, the probability of having k — 1 females in a 
single time step is Qk,k-i = k/2N, since one of the females has to be selected (probability 
k/N) and be replaced by a male offspring (probability 1/2). Similarly, the probability of 
changing to a state with k + 1 females is — (N — k)/2N and that of remaining with 

k females is = 1/2. The dynamics is similar to a random walk in the space of integers 
< k < N biased towards k = N/2: the closer A; is to the smaller the probability of 
k — > k — 1; the closer k is to iV the smaller the probability of k — > k + 1. This tends to 
stabilize the population, leading to small extinction probabilities. We define a generation 
by N time steps, which is the average time it takes for the entire population to be replaced. 

The dynamics of Pt{k) is governed by the following equations: 

P t+ i(0) = Pt(0) + ^P t (l) 

^+1(1) = ^Pt(l) + ^P(2) 

/V — k 4- 1 N k 4- 1 

Pt + i(k) = 2N + fi(fc-l) + w P«(AO + ^P«(fc + l), k = 2,3,...,N-2 
Pm(iV-l) = ^Pt(N - 2) + ^P t (N - 1) 



Pt+i(N) =^P t (N-l) + P t (N). 

(1) 

The probabilities Pt{k) define a vector of N + 1 components P t , in terms of which the 



master equation above becomes 

Pt+i = UP t . (2) 

The evolution matrix U and the transition matrix Q = U T are tridiagonal. This is a linear 
system that can be completely solved in terms of the eigenvalues and eigenvectors of U . 
However, because U is not symmetric, both right a r and left b r eigenvectors are needed. 
Moreover, since £V C/y = Yli^j,i = 1, U is a, stochastic matrix, having real eigenvalues 
satisfying Aj < 1 and 

JV 

^a r .? = l (3) 

r=0 

where the superscript T stands for transposition, the lower dot represents the diadic product 
and the normalization is set by of ■ 3j = 5ij. Using this property, the transition probability 
between an initial state with ko females and a state with k females after the time t can be 
written as 

N 

P(k, t; fc„, 0) = b r k a rk K (4) 

r=0 

where a r i is the i-th component of the r-th right eigenvector and similarly for b r j. Equation 
rtl]) is an example of a Markov chain, a discrete dynamical system where the transition 
probability between any two states depends only on the two states involved and not on the 
past history of the system. Markov chains are ubiquitous in genetics [TTH2"T] . but not so 
common in population dynamics. 

The dynamical system described by Eq. is related to the one dimensional motion of a 
Brownian particle subjected to an external force F(x), studied by Smoluchowski and Kac 
[28] . In this case the probability P(x, t) of finding the particle at position x at time t satisfies 
the diffusion equation 

dP d 2 P Id 

where D and / are the diffusion and friction coefficients. In the limit of large N we may set 
e = 1/N, Xk = ke and the time step to 5 and transform equation (JlJ) into a similar equation, 

OP e 2 d 2 P e d r , , ln . nl 

The 'diffusion coefficient' is e 2 /4S and the 'force' is harmonic, F(x) = — {x — 1/2), tending to 
restore the population towards x — 1/2 (or k = N/2). The main difference between Eqs.Q 
and ^ is that in the former the end points of the diffusion interval are absorbing states, 
lending the stationary solution of this equation unstable for any finite N. 



Since each of the N individuals can be either male or female and only one individual is 
replaced at each step, the dynamics described by Eqs.Q can also be mapped into a random 
walk to nearest neighbors in a hypercube in N dimensions This problem, it turn, is also 
related to the Ehrenfest model, where N numbered balls are placed into two boxes and at 
each time step one ball is chosen at random and moved to the other box [28] . Starting with 
all balls in one of the boxes and letting the system evolve corresponds to the classic model 
of a gas confined in one of two chambers and removing the wall separating the chambers. 
The state corresponding to all individuals of the same sex maps into two opposing corners 
of the hypercube and to all balls in the same box. 

In this paper tackle the problem of sex ratio fluctuations in population dynamics and 
make contact with these classic statistical models. Some of the known analytical results 
available in the literature cited above will be connected to the present calculations below. 



III. ABSORBING AND TRANSIENT STATES 



The eigenvectors corresponding to A = 1 completely determine the asymptotic behavior 
of the system, since the contributions of all the others die out at long times. 
The evolution matrix U is given explicitly by 



U 




There are two eigenvalues 1, with eigenvectors 

/l/2\ 



«o 








VV2/ 



and 



I V2 \ 





(7) 



(8) 



such that a + a at corresponds to P(0) and a — a at to P(N), both leading to extinction. 
The corresponding left eigenvectors are 



60 



/l\ 

1 
1 

1 

V 1 / 



and 



/ 1 \ 

bN,i 

bN,N-l 

V -1 / 



(9) 



The vector bjy does not have a simple form for finite N but the coefficients bj<r,k g° to zero 
for large N and k — 1,2 . . . , N — 1. Notice that this choice of vectors agrees with bj ■ Sj = 5ij 
in this sub-space. Using these vectors, equation Q can be re-written, for k 7^ 0, N, as 

N-l 



1 1 N ^ 

P(k, t; k , 0) = -S kfi + -S k ,N + 22 brk 



(10) 



r=l 



making explicit that extinction is the only asymptotic possibility, independent of the initial 

state k . 

The other eigenvalues and eigenvectors can be calculated from the non-trivial (N — 1) x 
(N — 1) part of U, delimited by the lines in equation ([?]), 



V = - + —W 
2 2N 



(11) 



where 



/ 
N-l 



W 



2 





3 



iV-2 4 















\ 






(12) 



40JV-2 

03 N-l 

\ 00 2 J 

If Aj are the eigenvalues of U, and /i, the eigenvalues of W, then 

1 



A, 



2 \ N 



(13) 
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IV. EIGENVALUES AND EIGENVECTORS FOR LARGE N 



Although W has a very simple structure, its eigenvalues cannot be calculated analytically 
for arbitrary N. It can be checked that the resulting polynomial for the eigenvalues \i gets 
more and more complicated as N increases. However, in the limit N — > oo, we find 



fH = N-2(i-l) 



1,2, 



(14) 



The largest eigenvalue, /ii = N, yields a new asymptotically stable state, with Ai = 1 and 
eigenvector 



a lk = ci exp 



-2N 



k 
N 



(15) 



corresponding to a symmetric distribution of males and females centered at k — N/2. In 
the limit of large N we may set e = l/N, x k = ke and obtain 

2" 



a(x) 




7Tf 



exp 



x 



(16) 



which is the stationary solution of equation ([6]). For finite N, however, Ai < 1 and the 
population inevitably goes extinct towards k = or k = N. Therefore, it is important 
to estimate how Ai tends to 1 as iV goes to infinity, a calculation presented in the next 
section. This provides an estimate of the time to extinction due to fluctuations in the sex 
determination. 

All other values of \ii lead to \ < Ai and do not contribute to the asymptotic state of 
the population. The eigenvectors are given by 



a ik = Ci exp 



-2N 



k 
N 



H^VWik/N- 1/2)), 



(17) 



where H m (x) are the Hermite polynomials and q are normalization constants. The appear- 
ance of Hermite polynomials in the solution can be traced to the diffusion equation ^ and 
its interpretation via Smoluchowski's Brownian motion. These results are demonstrated in 
the appendix [Aj 

V. LARGEST EIGENVALUE OF W FOR FINITE N 



Writing ^ = N — a, the eigenvalue equation det [W — /i]l] =0 becomes 

det [C + al] = det [C] det [1 + aC' 1 } = 0, 



(18) 



where C = W — Nl. Using 



det [1 + C- l a] = 1 + aTrlC- 1 ] + 0{a 2 ) (19) 
we obtain a = — {Tr[C -1 ]} , 

^ = N + {Trie- 1 ]}' 1 (20) 

and 



Ai = 1 + ^ {TriC- 1 ]}- 1 . (21) 



We show in Appendix [B] that 



JV-l 



TriC- 1 ] = V h- 1 ?"-"- 1 (22) 
where 

A = -Nf k . x -k(N-k+ l)/ fc _ 2 (23) 
with /o = 1 and /i = — iV. It turns out that, for large N, 

fk - 1 J N - k - 1 = -±-B(N, k)[l + 0(1/N)} (24) 

JN-l 

where B(N, k) is the binomial coefficient. In this approximation the sum can be easily 
performed and the result is 

A x = 1 - 2-<* +1 > w 1 - 2"^, (25) 

which is very accurate even for small N. Setting Ax = e _1 / Te we obtain the time to extinction 
as 2 N time steps. In terms of number of generations 

r e = 2 N /N. (26) 



The last approximation in Eq. (25), where we multiply by a factor 2, is justified because 
Eq.(24) is accurate only up to the scaling behavior of order 1/N. The factor 2 is obtained 
by fitting the numerical simulations (see figure [I]). The time to extinction can be related 
to Ehrenfest model. In this problem iV balls are placed in two boxes and a random ball is 
moved from its box to the other at each time step. Starting from an arbitrary state, the 
time it takes for finding all the balls in one box is of order 2 N [28]. This is also the Poincare 
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recurrence time, which is the time it takes to return to the state with all balls in the same 
box, having started there. 

Considering the meta-stable probability distribution P(k) = 2~ N B(N, K) corresponding 
to Ai, one could ask how long it takes for an arbitrary initial state to reach P(k). This is 
given by the next eigenvalue, A 2 = 1 — 1/N, with the associated relaxation time of a single 
generation: 

r r = 1. (27) 

This is the second important time scale of the problem, much shorter than r e . This result 
is analogous to that obtained in [29] for the time taken by a particle to reach the stationary 
distribution on the hypercube under a nearest neighbor random walk, which is of order 
iVlogiV. 

As a simple example of these time scales, a population with N = 20 individuals starting 
with k — 10 females fluctuates according to the binomial distribution after 1 generation, or 
20 time steps, but goes extinct only after about r e = 2 20 /20 ~ fifty thousand generations 
(or a million time steps). Figure [I] shows a comparison between numerical calculations and 
the theoretical prediction of r e , showing good agreement even for small values of N. 

VI. ROBUSTNESS AGAINST CHANGES IN OFFSPRING SEX RATIO 

In some species the birth of one of the sexes is favored over the other. A well studied 
example is the infection caused by the bacteria Wolbachia, that kills males in some arthropod 
species [30]. Sex determination can also be influenced by temperature in several species of 
reptiles, like the snow skink lizard |31j . Changes in the sex ratio of offspring not only shifts 
the distribution of males and females in the population but also affects the time to extinction 
due to random fluctuations. In this section we discuss the time to extinction assuming that 
the probability of a female offspring is p = 1/2 + s and that of a male is 1 — p = 1/2 — s. 

As we demonstrate below, the inclusion of bias in the sex ratio at birth complicates 
the dynamics and only approximate solutions for the time to extinction can be derived. 
However, the problem can be easily solved for the Wright-Fisher model, where generations 
are non-overlapping and constructed from the previous one by independent random choice 
of males and females. The Wright-Fisher model describes, for instance, annual plants, where 
the entire population dies in the winter and its replaced anew in the spring. The Moran 
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FIG. 1. (Color online) Time to extinction as a function of N. The black points show the average 
and mean square deviation obtained with 100 replicates of simulations. The red line shows the 



theoretical result, Eq. (26) 



model, on the other hand, is appropriate for perennial plants. The probability of k females 
is given by 

(l/2 + s) k (l/2-s) N ~ k B(N,k) (28) 

and the probability of extinction is, therefore (1/2 + s) N + (1/2 — s) N , which reduces to 
2-(iV-i) f or s _ q fj^g time t extinction is the inverse of this probability and is already 
given in terms of number of generations: 2^ Ar_1 - ) = (2 N /N) x (N/2). Although this can be 
taken as a first estimate for the time to extinction in the Moran model, it overestimates it 



by a factor (N/2) (see Eq.(26)). The factor N comes from updating the entire population 



at once. Numerical simulations indicate that the factor 1/2, valid for s=0 only, becomes 
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smaller for more extremes values of s. 

In the case of the Moral model, the transition matrix elements are generalized to 



n k>k = ^(1 + 2s) + ^(1 - 2s) = ±(N + 2s{2k - N)) 



(29) 



N-k, 



2N (1 + 2S ) 

and a master equation similar to equation (JTJ can be written. In the extreme cases s = ±1/2 
all the eigenvalues of the evolution matrix U can be calculated analytically and the largest 
non-unit eigenvalue is Ai = 1 — l/N, so that the time to extinction in approximately N, i.e., 
1 generation. 

In the limit of large N, the master equation can also be written as a diffusion equation 



dP _ d 
dt dx 



D(x, s 



dP 

dx 



ed_ 
5 dx 



[F(x,s)P] 



(30) 



with a 'space-dependent' diffusion coefficient D(x, s) = (e 2 /45)(l + 2s — 4sx) and F(x, s) 
— (x — 1/2 — s — es). The stationary solution is 



P(x) = Aexp 



4sx + [1 - 4(1 + e)s 2 ] ln(l + 2s - 4sx) 
4es 2 



(31) 



where A is a normalization constant. Notice the symmetry of these equations with respect 
to the change s — > —s and x — > 1 — x. 

Differential equations for the continuous limit of the coefficients a k and b k can also be 
obtained and the corresponding solutions are \ii — N — 2(i — 1) with 



bi(x) = CiHi-x 



2N 



I -As 2 



(x-l/2-s) 



and 



2N(x - 1/2 - s) 2 . , , , 
ai(x) = d^exp <^ — >bi{x), 



(32) 



(33) 



1 -4s 2 

where Hi(x) are the Hermite polynomials. Although the solution for a\(x) looks rather 



different from (31), they are very similar for large N. Equation (33) for i = 1, corresponding 
to [Ax = N, is just a Gaussian centered at x = 1/2 + s with variance 1 — 4s 2 , as it should be. 
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FIG. 2. (Color online) Time to extinction as a function of s for a population with N = 20 starting 
with 10 males and 10 females. The line with square symbols shows the result of simulations 
obtained from 10,000 realizations for each value of s. The thick red line shows the approximation 
r 2 = a2 N «ff/N eff for a = 2. 



The correction to the \i\ is still given by equation (20) but now with 



N-1 , ( , , N-1 



k=i 



k=i 



where 



f k {s) = -[N + 2s(N - 2k)]f k „ 1 {s) - k(N — k + 1)(1 - 4 S 2 )/ fc _ 2 (s) 



(34) 



(35) 



with /o = 1 and f\ = —N + 2s(N — 2). Explicit evaluation of the coefficients d(N,k,s), 
however, is much more complicated and an analytic expression is not available. Numerical 
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simulations show that, for large N, we may approximate 
where 

A(N, s) = 2 Af ( 1 - L4 W) 2 (37) 
was obtained by fitting the amplitude of the coefficients. Further approximating the sum 



over k in equation (34) by an integral we obtain 



-/vr-'| i \ - \r°J° 2) (38) 



so that the time to extinction becomes 



Te = 2y ?r(1 2iV 4g2) 2 jV(1 - L4|s|)2 ~ 2""' (39) 

with N e ff = N(l — 1.4|s|) 2 . Figure 2 shows the time to extinction as a function of s 
for iV = 20. The solid line shows the result of simulations and the red line a fit with 
T2 = a2 Neff /N e ff. The goodness of the fit suggests that the expression is correct and that 
the main effect of s is to change the population size to an effective value. 

Our calculations have shown that producing males and females offspring at the same 
proportion is a good strategy even for very small populations. If, however, the birth of one 
of the sexes is favored, not only the balance between males and females is altered but the time 
to extinction might decrease dramatically. For small deviations, however, the exponential 
character of the extinction time remains. 
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Appendix A: Eigenvalues and eigenvectors of W for large iV 

Setting ciio = = ho = biN = we obtain the following recurrence relations for the i-th 
right and left eigenvectors of W: 

(k + l)a ik+ i + (N - k + l)a iJfc _i - = (Al) 
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and 

(N - k)b ik+1 + kb ik _i - Hib ik = 0. 
It can be checked that these components are related by 



a ik = 2~ N Bi(N, k)b 



ik ■ 



(A2) 



(A3) 



where Bi(N,K) is the binomial coefficient. 



Dividing equation (A2) by A" and defining the function bi(x) such that bi(x k ) = bi k leads 



to 



(1 - x)bi(x + e) + xbi(x - e) - ^bi(x) = 0, 



(A4) 



where e = 1/N. 

For N ^> 1, we can approximate b~i{x ± e) = b~i(x) ± tb'^x) + ^-b"(x). Neglecting terms of 
0(e 3 ), equation (A4) is turned into a differential equations for b(x): 



bf- + b' l e(l-2x) + b l (l-|) =0. 



(A5) 



Taking the limit A^ — > oo (e — > 0) we see that — N is an eigenvalue. The corresponding 
left and right eigenvectors are given by 



>ik 



-N 



for k = 1---N- 1. 



A;!(A^ — £;)! 

To obtain the remaining eigenvalues and eigenvectors for large A^ we define t 



1/2) and the function g{t) = b (^= ). Accordingly, equation (|A5|) becomes 



g" - 2td + 2d i9i = 



(A6) 
(A7) 



where di = (N — yUi)/2. Equation (A7) is satisfied by the Hermite polynomials if di 



were integers. To see that this is indeed the case, note that the boundary conditions are 
6(0) = 6(1) = 0, where the argument is x = t/y/2N + 1/2. Therefore, the corresponding 
boundary conditions for g(t) are g(±^N/2) = or, for A^ — > oo, g(±oo) = 0. In order 

to prevent g(t) from diverging we must set di — 0, 1, 2, The procedure is similar to the 

quantization of the harmonic oscillator in quantum mechanics. In this limit the eigenvalues 
and eigenfunctions are given by 

fH = N-2(i- 1) (A8) 
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and 

gi (t) = CiHi-^t) (A9) 

or 

b ik = c i H i - 1 (V2N(k/N-l/2)) (A10) 

where % — 1, 2, 3, . . . and the q are normalization constants. 

For the right-eigenvectors we can approximate the binomial by the normal distribution 
Bi(p, N) PS N(fi, a) by taking p = 1/2, \x = Np and cr = Np{l — p). We obtain 

a ik = c'texpftNik/N - l/2) 2 ]H i _ 1 (V2N(k/N - 1/2)). (All) 



Appendix B: The trace of C 



-i 



The matrix C is given by 

( -N 2 
N — 1 —N 3 



C = 



V 



o 
o 






N — 2 —N 4 















4 -N N - 2 
3 












-N N — 1 
2 -N J 



(Bl) 



and we only need the diagonal elements of C 1 , which can be obtained by Laplace's formula: 

Cofactor(k, k) 



[C 



-ii 



kk — 



(B2) 



Det[C] 

The Cofactor(ij) is the determinant of the auxiliary matrix obtained by removing the j-th 
line and i-th row of C, multiplied by (— 

The removal of the k-th line and row of C divides the remaining matrix into two decoupled 
blocks. In order to deal with these blocks we recursively define the matrix C^-k to be 
the matrix C^-k+i with the first line and first row removed. This definition holds for 
k — 2, . . . , N — 1, with the initial condition Cn-i = C. Calling fN-k = Det[CN-k] it is easy 
to see that 

Cofactor(k,k) = fk-ifw-k-i- (B3) 
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Applying the Laplace rule to the determinant f k it can be checked that it satisfies the 
recurrence relation 

/* = -Nf k . x — k(N — k + l)/ fc _ 2 (B4) 

with /o = 1 and f\ = — AT. 



In order to calculate (B3) it is useful to define = (— 1) f k -i/k\. In terms of g k the 



recurrence relation becomes 



(k + l)g k+1 - Ng k + (N -k + l)g k ^ = 



(B5) 



with g = and g\ = 1. We obtain 



fk-lfN-k-l 9k 9N-k I 



lA'/i- 



AT-1 



^ 5(iV,fc) 



(B6) 



where B(N,k) is the binomial coefficient. Since (B5) is the relation satisfied by B(N,k) 
itself, it is reasonable to assume that 



g k = ^B(N,k) 



(B7) 



where the factor 1/N guarantees the initial condition g\ — 1. This, however, is only an 
approximation, since it gives go = 1/N and not zero. However, for large ./V it suffices for 



obtaining the first order correction to the eigenvalue. Replacing (B7) into (B6) we obtain 

B(N, k) 



\c-\ 



kk 



N 



(B8) 
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